## File Description: Examine wells from old gas states wells and compare to ONRR data
# Author: Brian Prest

#****# Read Data, Compute Dates and Well Types -----------------------------------------------------------------

# Read in headers. Data is split into 4 files due to download restrictions.
dt.oldgas = fread(di.data.folder%&%'/OLDGAS-2020-04-14--Top 4 Gas States Pre 1990/HPDIHeader.csv', drop=221, na.strings="(N/A)") # data tables thinks there's a 221th col?
names(dt.oldgas) <- tolower(names(read.csv(di.data.folder%&%'/OLDGAS-2020-04-14--Top 4 Gas States Pre 1990/HPDIHeader.csv', nrows=1))) # data tables failed to read colnames
dt.oldgas

dt.oldgas[, first_prod_date := as.Date(first_prod_date, format='%m-%d-%Y')]
dt.oldgas[, first_prod_month := first_prod_date - day(first_prod_date) + 1]

## Label wells as oil or gas
dt.oldgas[prod_type=='OIL (CYCLIC STEAM)', prod_type := 'OIL']
dt.oldgas[!(prod_type %in% c('OIL','GAS','CBM')), prod_type := 'Other']
dt.oldgas[prod_type=='CBM', prod_type:='GAS']
dt.oldgas[prod_type=='OIL', prod_type:='Oil']
dt.oldgas[prod_type=='GAS', prod_type:='Gas']

dt.oldgas[prod_type=='Other' & gor_cum>oil.gor.p90, prod_type:='Gas'] # Classify these as Gas
dt.oldgas[prod_type=='Other' & gor_cum>0, prod_type:='Oil'] # For the rest, label them as oil

# Classify the ND Other wells as oil
dt.oldgas[prod_type=='Other' & state=='ND', prod_type:='Oil']
# Classify the remaining ones as gas
dt.oldgas[prod_type=='Other', prod_type:='Gas']

# Federal lands overlay
locs = dt.oldgas[, .(longitude, latitude)]


coordinates(locs) <- c('longitude','latitude')
class(locs)
class(vessel)

proj4string(vessel)

proj4string(locs) <- proj4string(vessel)

over.out = over(locs, vessel)
over.out = as.data.table(over.out)
names(over.out) = tolower(names(over.out))
over.out = cbind(as.data.table(locs), over.out)

dt.oldgas[, offshore:=0] # these states are landlocked 

dt.oldgas = cbind(dt.oldgas, over.out)
dt.oldgas[, fed.land := 1*!is.na(admin1)]


dt.oldgas.panel = fread(di.data.folder%&%'/OLDGAS-2020-04-14--Top 4 Gas States Pre 1990/HPDIProduction.csv')
names(dt.oldgas.panel) <- tolower(names(read.csv(di.data.folder%&%'/OLDGAS-2020-04-14--Top 4 Gas States Pre 1990/HPDIProduction.csv', nrows=1))) # data tables failed to read colnames
dt.oldgas.panel[, days:=NULL]
dt.oldgas.panel[, wcnt:=NULL]
dt.oldgas.panel[, wtr:=NULL]
dt.oldgas.panel[, prod_date := as.Date(prod_date, format='%m-%d-%Y')]

dt.oldgas.panel = merge(dt.oldgas[,.(entity_id, first_prod_month, state, prod_type, fed.land, county)], dt.oldgas.panel, by='entity_id', all.y=FALSE)

# 1250.5+565.7 = 1816.2 bcf of gas produced in 2018, 1250 of which is on fed land.
# This still understates ONRR's fed production of 1360 somewhat, but it overstates total EIA production of 1720
dt.oldgas.panel[year(prod_date)==2018, sum(gas)/1e6, by=.(fed.land)]
dt.oldgas.panel[year(prod_date)==2018, sum(gas)/1e6, by=.(fed.land, new=year(first_prod_month)>=1990)]
# Note that 259+104=364 out of the 1816 (20%) of WY gas production is from wells predating
# 1990. So the 1990 cut understates things a bit.
# It must be that the gas stripper wells.

dt.oldgas.panel[, offshore:=0] # these states are all landlocked (WY,CO,UT,NM)

save('dt.oldgas.panel', file=working.data.dir%&%'merged_dt_panel_oldgas_'%&%today()%&%'.RData')
oldgas_loc = working.data.dir%&%'merged_dt_panel_oldgas_'%&%today()%&%'.RData'

tot.prod.by.fed.old.gas.states = dt.oldgas.panel[, .(liq=sum(liq)/(365.25/12)/1e6, gas=sum(gas)/(365.25/12)/1e6), 
                               by=.(prod_date, fed.land, offshore)]
tot.prod.by.fed.old.gas.states = tot.prod.by.fed.old.gas.states[order(fed.land, offshore, prod_date)]
tot.prod.by.fed.old.gas.states[, land.type := ifelse(fed.land==1, 'Federal','NonFederal')]
tot.prod.by.fed.old.gas.states[, offshore := ifelse(offshore==1, 'Offshore','Onshore')]
tot.prod.by.fed.old.gas.states = dcast(tot.prod.by.fed.old.gas.states, prod_date ~ land.type + offshore, value.var=c('liq','gas'))
tot.prod.by.fed.old.gas.states[, liq_NonFederal_Offshore:=0]
tot.prod.by.fed.old.gas.states[, gas_NonFederal_Offshore:=0]
tot.prod.by.fed.old.gas.states[, liq_Federal_Offshore:=0]
tot.prod.by.fed.old.gas.states[, gas_Federal_Offshore:=0]
tot.prod.by.fed.old.gas.states[, liq_tot := liq_Federal_Onshore + liq_Federal_Offshore + liq_NonFederal_Onshore + liq_NonFederal_Offshore]
tot.prod.by.fed.old.gas.states[, gas_tot := gas_Federal_Onshore + gas_Federal_Offshore + gas_NonFederal_Onshore + gas_NonFederal_Offshore]

save('tot.prod.by.fed.old.gas.states', file=working.data.dir%&%'/old_gas_states_total_'%&%today()%&%'.RData')

rm(dt.oldgas)
rm(dt.oldgas.panel)
